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Abstract 

We present two algorithms for computing the geodesic distance between phylogenetic trees in 
tree space, as introduced by Billera, Holmes, and Vogtmann (2001). We show that the possible 
combinatorial types of shortest paths between two trees can be compactly represented by a 
partially ordered set. We calculate the shortest distance along each candidate path by converting 
the problem into one of finding the shortest path through a certain region of Euclidean space. 
In particular, we show there is a linear time algorithm for finding the shortest path between 
a point in the all positive orthant and a point in the all negative orthant of M'^ contained in 
the subspace of consisting of all orthants with the first i coordinates non-positive and the 
remaining coordinates non-negative for < i < fc. The resulting algorithms for computing the 
geodesic distance appear to be the best available to date. 

1 Introduction 

Phylogenetic trees, or phylogenies, are used throughout biology to understand the evolutionary 
history of organisms ranging from primates to the HIV virus. Outside of biology, they are used in 
studying the evolution of languages and culture, for example. Often, reconstruction methods give 
multiple plausible phylogenetic trees on the same set of taxa, which we wish to compare using a 
quantitative distance meaure. We would also like to have a statistical framework to better evaluate 
the generated trees. The tree space of Billera, Holmes, and Vogtmann [T and its correspond- 
ing geodesic distance measure address both of these issues. In this paper, we give two practical 
algorithms for computing this distance. 

There are many different algorithms to construct phylogenetic trees from biological data ( |10j 
and its references), but their accuracy can be affected by such factors as the underlying tree shape 
[13] or the rate of mutation in the DNA sequences used [13] • To compare these methods through 
simulation, or to find the likelihood that a certain tree is generated from the data, researchers need to 
be able to compute a biologically meaningful distance between trees [11] . Several different distances 
between phylogenetic trees have been proposed (e.g. [7], [11], [12], [21], [22]). For example, the 
Nearest Neighbor Interchange (NNI) distance ([2T] and j29j). which counts the number of rotations 
between trees, is considered one of the best distance measures, but is NP-hard to compute |8j, 
making it impractical. 

In response to the need for a distance measure between phylogenetic trees that naturally in- 
corporates both the tree topology and the lengths of the edges, Billera et al. [3] introduced the 
geodesic distance. This distance measure is derived from the tree space, %i, which contains all 
phylogenetic trees with n leaves. The tree space is formed from a set of Euclidean regions, called 
orthants, one for each topologically different tree. Two regions are connected if their corresponding 
trees are considered to be neighbours. Each phylogenetic tree with n leaves is represented as a 
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point within this space. There is a unique shortest path, called the geodesic, between each pair of 
trees. The length of this path is our distance metric. 

The most closely related work is by Staple j26j and Kupczok et al. fTS], who developed al- 
gorithms to compute the geodesic distance based on the notes of Vogtmann [28]. Both of these 
algorithms are exponential in the number of different edges in the two trees. Although Kupczok et 
al. developed their algorithm GeoMeTree independently, it can be considered a direct improve- 
ment to the algorithm of Staple. We show in Section 5 that our algorithm performs significantly 
better than GeoMeTree, although it is still exponential. A polynomial time, \/2-approximation 
of the geodesic distance was given by Amenta et al. [T] . 

Our primary contribution is two algorithms for computing the geodesic distance between two 
phylogenetic trees. We show these algorithms are significantly faster than the only explicit algorithm 
published to date. Furthermore, two main ideas were developed to construct these algorithms. First, 
the candidate shortest paths between trees can be represented as an easily constructible partially 
ordered set, giving information about the combinatorics of the tree space. Second, we can find the 
length of each candidate shortest path by translating the problem into one of finding the shortest 
path through a region of a lower dimensional Euclidean space. The solution to this new problem is 
a linear algorithm for a special case of the shortest Euclidean path problem in M" with obstacles. 
Since the general problem is NP-hard for dimensions greater than 2, this result is also of interest 
to computational geometers. These two ideas can be combined using dynamic programming or 
divide and conquer methods to significantly reduce the search space, and thus make this distance 
computation practical for some biological data sets of interest. 

The remainder of this paper is organized as follows. In Section 2, we describe the tree space and 
the geodesic distance between phylogenetic trees. The problem of finding the geodesic distance has 
both a combinatorial component, which is investigated in Section 3, and a geometric component, 
which is covered in Section 4. More specifically, we introduce a combinatorial framework in Section 
3, which represents the candidate shortest paths between trees by an easily constructible partially 
ordered set (Theorem 3.7). In Section 4, we translate the problem of calculating the length of a 
candidate shortest path into a problem in Euclidean space (Theorem 4.4), and then show that this 
Euclidean problem can be solved in linear time. Section 5 combines the ideas of Sections 3 and 4 
using dynamic programming or divide and conquer techniques to present two complete algorithms. 



2 Tree Space and Geodesic Distance 

This section describes the space of phylogenetic trees, 7^, and the geodesic distance. For further 
details, see [3]. A phylogenetic tree, or just tree, is a rooted tree, whose leaves are in bijection with a 
set of labels X representing different organisms. For this paper, let X = {1, n}. We often treat 
the root as a leaf, called 0. We consider both bifurcating (or binary) trees, in which each interior 
vertex has degree 3, and multifurcating (or degenerate) trees, in which this is not the case. 

A split A\B IS a partition of X U {0} into two non-empty sets A and B, where X is the leaf-set 
of some tree T and is its root. A split is in T if it corresponds to some edge in T, in that one 
block of the partition consists of all the leaves below that edge, while the other block consists of 
the remaining leaves and the root. We say this split is induced by that edge in T. For example, in 
Figure [T] the split induced by the edge 63 partitions the leaves into the sets {2,3} and {0, 1,4,5}. 
We will refer to splits induced by an edge ending in a leaf as trivial splits, and to all other splits as 
simply splits. A split of type n is a partition of the set {0, 1, n} into two blocks, each containing 
at least two elements. Let Et be the set of (non-trivial) splits of tree T. A £ Et is a set of splits 
in T, then let T/A be the tree T with each edge that induces an element of A contracted. 
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Figure 1: The spht induced by the edge 63. 



Two sphts e = X\X' s.nAe' = Y\Y' are compatible if one of X n X n F', X' n F and X' n Y' 
is empty. Equivalently, two sphts are compatible if their inducing edges can exist in the same 
phylogenetic tree. For example, in Figure [T] the split 63 = {2, 3}|{0, 1, 4, 5} is compatible with 
the split 62 = {2, 3, 4}|{0, 1, 5}, because {2,3} n {0,1,5} = 0. However, 63 is incompatible with 
/ = {1, 2}|{0, 3, 4, 5}. Two sets of mutually compatible splits of type n, A and B, are compatible if 
^ U i? is a set of mutually compatible splits. 

Each edge, and hence split, e E Et-^ is also associated with a non-negative length |e|r- For 
example, this length could represent the number of DNA mutations that occurred between specia- 
tion events. Two splits are considered the same if they have identical partitions, regardless of their 

lengths. For any A C Et, let ||^|| = w EeeA l^lr- 



2.1 Tree Space 

We now describe the space of phylogenetic trees, 7^, as constructed by Billera et al. [3]. It is 
homeomorphic, but not isometric, to the tropical Grassmannian [21] and the Bergman fan of the 
graphic matroid of the complete graph [2j. This space contains all bifurcating and multifurcating 
phylogenetic trees with n leaves. In this space, each tree topology with n leaves is associated with 
a Euclidean region, called an orthant. The points in the orthant represent trees with the same 
topology, but different edge lengths. These orthants are attached, or glued together, to form the 
tree space. 

We do not use the lengths of the edges ending in leaves in the definition of tree space, but can 
easily include them by considering geodesies through 7^ x M" , as noted in Billera et al. [3]. 

Any set of n — 2 compatible splits corresponds to a unique rooted phylogenetic tree topology 
( |23| Theorem 3.1.4]). For any such split set Et corresponding to tree T, associate each split with 
a vector such that the n — 2 vectors are mutually orthogonal. The cone formed by these vectors is 
the orthant associated with the topology of T. Recall that the fc-dimensional (nonnegative) orthant 
is the non-negative part of M'^, denoted M^p. A point (xi, Xn_2) in represents the tree 

containing the edge associated with the z-axis that has length Xi, for all 1 < i < n — 2, as illustrated 
in Figure p (a) [ If Xi = 0, then the tree is on a face of the orthant. In this case, we will say the tree 
does not contain the edge associated with the i-axis. The trees on the faces of each orthant have 
at least one edge of length 0. Furthermore, two orthants can share the same boundary face, and 



thus are attached. For example, in Figure p( a) the trees Ti and T[ are represented as two distinct 
points in the same orthant, because they have the same topology, but different edge lengths. The 
tree Tq has only one edge, ei, and thus is a point on the ei axis. 



Notice that although Figure 2(a) is drawn in the plane, it actually sits in M^, with each of the 
axes or splits corresponding to a different dimension. In general, 7^, sits in M-^, where M is the 
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number of possible splits. However, as no point in 7^ has a negative coordinate in M-^ , we will 
often let the positive and negative part of an axis correspond to different splits. 



geodesic 





— — = geodesic 



(a) Two orthants in 74. 



(b) Both edge length and tree topology determine the 
geodesic. 



Figure 2: The geometry of tree space. 

For any set A of compatible splits with lengths, let T(A) represent the tree containing exactly 
the edges that induce the splits in A and all trivial splits. The lengths of the edges in T{A) 
correspond to their respective lengths in A. Let 0{A) be the orthant of lowest dimension containing 
T{A). For any t > 0, let t ■ A he the set of splits A whose lengths have all been multiplied by 
t. If A and B are two sets of mutually compatible splits of type n, such that AD B is also a 
set of mutually compatible splits, then we define the binary operator + on the orthants of 7^ by 
0{A) + 0{B) = 0{A U B). For any union of disjoint orthants, uf^oC'(^i), where 5 is a set of 
mutually compatible splits of type n such that B and Ai are compatible sets for all < i < /c, 
define {u^^QO{Ai)) + 0{B) = U^Lq {0{A,) + 0{B)) = U^^f)0{Ai U B). UAnB = 0, we also use 
the direct sum notation ©. 



2.2 Geodesic Distance 

There is a natural metric on 7^. The distance between two trees in the same orthant is the Euclidean 
distance between them. The distance between two trees in different orthants is the length of the 
shortest path between them, where the length of a path is the sum of the Euclidean lengths of the 
intersections of this path with each orthant. For any trees Ti and T2 in 7^, the geodesic distance, 
d{Ti,T2), between Ti and T2 is the length of the geodesic, or locally shortest path, between Ti and 
T2 in Tn. Billera et al. defined this distance, and proved that 7^ is a CAT(O) space p] Lemma 4.1], 
or has non-positive curvature [5], and thus the geodesic between any two trees in 7^ is unique. 

For example, in Figure 2(a) the geodesic between the trees Ti and T2 is represented by the 
dashed line. Figure 2(b) depicts 5 of the 15 orthants in 74. This figure also illustrates that the 
edge lengths, in addition to the tree topologies, determine through which intermediate orthants the 
geodesic will pass. 



2.3 The Essential Problem 

The problem of finding the geodesic between two arbitrary trees in 7^ can be reduced in polynomial 
time to the problem of finding the geodesic between two trees with no splits in common. This is 
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the problem considered in Sections 3 and 4. Furthermore, the lengths of the pendant edges can 
easily be included in the distance calculation, if desired. 

Vogtmann [28^ proved the following theorem, which explains how to decompose the problem of 
finding the geodesic when the trees share a common split. An alternative proof is given in [19]. Let 



Ti and T2 be two trees with a common split e = X\Y, where G X, as shown in Figure 3(a) For 



i E {1,2}, let T/^ be the tree Tj with edge e and any edge below it in Ti contracted. That is, any 
edge a = X'\Y' such that X' C 1" or y' C y is contracted, as shown in Figure |3(b)| For i e {1,2}, 
let be the tree Tj formed by contracting edge e and all edges not contracted in T^. That is, 
any edge b = X'\Y' such that X' d X 01 Y' d X \s contracted, as in Figure [3(c) [ 








(a) Tree Ti. (b) Tree T/*. (c) Tree tP . 

Figure 3: Forming the trees and Tj^ from Tj for i G {1,2}. 

Theorem 2.1. If Ti and T2 have a common split e, and and are as described in the above 
paragraph for i G {1,2}, then d{Ti,T2) = ^d{T^,T^Y + d{Tf,T^)^ + {\e\T, - \e\T2f . 

As noted above, the length of the edges ending in leaves can be included in the distance calcula- 
tions by considering the product space 7^ x M", and the shortest distance, d;(Ti,T2), between the 
trees in this space. In this case, if the length of the edge to leaf i in tree T is \li\T for all 1 < i < n, 

then di{Ti,T2) = ^d{Ti,T2)^ + {\li\T, - Iklr^f- 

Therefore, the essential problem is as follows, and we devote the rest of this paper to it. 

Problem. Find the geodesic distance between Ti and T2, two trees in 7^ with no common splits. 



3 Combinatorics of Path Spaces 

The properties of the geodesic imply that it is restricted to certain orthants in the tree space. 
In this section, we model this section of tree space as a partially ordered set (poset), called the 
path poset, in which each element corresponds to one of the orthants in the tree space. This poset 
enables us to enumerate all orthant sequences that could contain the geodesic, as each such orthant 



sequence, called a path space, corresponds to a maximal chain in the path poset by Theorem 3.7 
For this section, assume that Ti and T2 are two trees in 7^ with no common splits. That is 

n = 0- 



3.1 The Incompatibihty and Path Partially Ordered Sets 

We now define the incompatibility poset, which depicts the incompatibilities between splits in Ti 
and T2. It will be used to construct the path poset. To define these posets, we introduce the 
following two split set definitions. 
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Let A and B be two sets of mutually compatible splits of type n, such that AnB = 0. Define 
the compatibility set of A in B, Cb{A), to be the set of splits in B which are compatible with all 
splits in A. Define the crossing set of A in B, Xb{A), to be the set containing exactly those splits 
in B which are incompatible with at least one split in A. We will use the set of mutually compatible 
splits and the tree containing exactly those splits interchangeably here, thus writing Ct{A) instead 
of CEj,iA) and Xt{A) instead of Xet,{A). 

If D is a set of mutually compatible splits of type n such that D <^ A, then: 

1. Cb{A) C Cb{D) {opposite monotonicity of the compatibility set), 

2. Xb{D) C Xb{A) {monotonicity of the crossing set), 

3. Cb{A) and Xb{A) partition B {partitioning). 

A preposet or quasi-ordered set is a set P and binary relation < that is reflexive and transitive. 
See \25\ Exercise 1] for more details. Define the incompatibility preposet, P{Ti,T2), to be the 
preposet containing the elements of Et2-, ordered by inclusion of their crossing sets. So, for any 
/, /' S Et2-, f < f P{Ti,T2) if and only if X'r^{f) C XT^{f'). Define the equivalence relation 
f ^ f \i and only f < f and /' < /. Thus, all the splits in an equivalence class have the same 
crossing set, which we define to be the crossing set of that equivalence class. 

Definition 3.1. The incompatibility poset, P{Ti,T2), consists of the equivalence classes defined by 
~ in the preposet P{Ti,T2) ordered by inclusion of their crossing sets. 

Generally, we will be informal, and treat the incompatibility poset as the elements of 
ordered by inclusion of their crossing sets. When we say two elements of P{Ti,T2) are equiva lent, 
we mean that formally they are in the same equivalence class in the preposet P(Ti, T2). Figure 4(c) 
gives an example of an incompatibility poset. 





2 5 6 



fa - {ei, 62, 63, 64} 
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_L 
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h - {63, 64} ^f4^ 
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U - {64} 







(a)Treeri. (b) Tree Ta. (c)P(Ti,r2) {A) K{T^,T2) 

Figure 4: Example of an incompatibility poset and a path poset. 
For any A G Et2 , define A £ Et2 by 

A^A = {feET2: XT,{f) C XtAA)}. 

Note that by definition, X'ri{A) = Xtj^{A). The map X 1-^ X is a closure operator on a set / if 
for every subset X C I, it is extensive {X C X), idempotent {X = X), and isotone (if X C Y, 
then X C Y) [4 . From the definitions and the monotonicity of crossing sets, A ^ A is a closure 
operator on Et2- 
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Definition 3.2. The path poset ofTi to T2, K{Ti,T2), is the closed sets of ordered by inclusion. 



The path poset represents the possible orthant sequences containing the geodesic between Ti 
and T2, and we next make clear this correspondence. The path poset is bounded below by 0, and 
above by E j-y. It is a sublattice of the lattice of order ideals of P{Ti,T2), but need not be graded 



pTQJ. Figure 4(d) gives an example of a path poset. For simplicity, we often just write /1/3 instead 



of {/i, Is}, for example. 
3.2 Path Spaces 

The geodesic is contained in some sequence of orthants connecting the orthants containing Ti and 
T2. Billera et al. |31 defined a set of orthant sequences, such that at least one of them contains the 
geodesic. We call such orthant sequences path spaces. We characterize all maximal path spaces 



in Theorem 3.4 and show that they are in one-to-one correspondence with the maximal chains in 



K(Ti,T2) in Theorem 3.7 



Definition 3.3. For trees Ti and T2 with no common splits, let Eq Z) Ei Z) ... D E^^i D E^^ and 
-Fq C -Fi C ... C Fk-i C Fk be sets of splits such that Ei and Fi are compatible for all < i < A:, 
Eq = Fk = Fx2i and E^ = Fq = 0. Then U^^QO{Ei U Fi) is a path space between Ti and T2. 

Note that the inclusions are strict in this definition. A path space is a subspace of 7^ consisting 
of the closed orthants corresponding to the trees with interior edges Ei U Fi for all < i < k. To 
simphfy notation, let d = 0{Ei U Fi) and 0[ = 0{E[ U Fl). The intersection of Oi and Oj+i is 
the orthant 0{Ei^i U Fi). If the i^^ step transforms the tree with splits Ei-i U Fi^i into the tree 
with splits Ei U Fi, then at this step we remove the splits Ei-i\Ei and add the splits 

A path space is maximal if it is not contained in any other path space. Since |i3, Proposition 
4.1] proves that the geodesic is contained in a path space, it must be contained in some maximal 
path space. We now characterize the maximal path spaces using split compatibility. 

Theorem 3.4. The maximal path spaces are exactly those path spaces U^^QOi such that: 

1. Ei = Cr^Fi), for allO<i<k. 

2. Fi = Cr^iEi), for allO<i<k. 

3. for all 1 < i < k, the splits in Fj\Fj_i are minimal and equivalent elements in the incompat- 
ibility poset P {T{Ei^i),T{ET2\Fi-i)) 

Notice that Conditions 2 and 3 imply that if / G Fi\Fi-i, then any other split /' equivalent 
with / is also in Fi\Fi^i, and these are precisely the splits in Fi\Fi^i. In contrast, an arbitrary 
path space has Ei C CxiiFi) and Fi C CT2{Ei), but not necessarily equality. 



Before giving a proof of Theorem 3.4 we define a relaxed path space. The only difference 
between a relaxed path space and a path sapce is that the inclusions need not be strict. We then 
show that any relaxed path space can be expressed as a path space. 

Definition 3.5. For trees Ti and T2 with no common splits, let Et-^ = Eq ^ Ei Z) ... D E^-i 5 
Ek = 0, and = Fq Q Fi (Z ... Q F^-i ^ Fk = Et2 be sets of splits such that Ei and Fi are 
compatible for all < i < /c. Then U^^qCj is a relaxed path space between Ti and T2. 

Lemma 3.6. Let S = U^^qCj he a relaxed path space. Then S is also a path space. 

Proof. IfS = uf=oCi 

is a relaxed path space, then one of the following cases holds: 
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Case 1: For some < j < k, Ej = -Ej+i and Fj C Fj 
Then Oj C Oj+i, and so S = (u{ZoOA U 



• Case 2: For some < j < k, Ej 5 -E'j+i and Fj = -Fj+i 



Then D O^+i, and so 5 = U^=o^0 U U^^ .+a©, 



Reindex the E'j's and Fj's in this new expression for 5 so that their indices are consecutive from 
to k — 1. Redefine k to be k — 1. While Case 1 or Case 2 holds, repeat this process. Since < A; < oo 
and we reduce k by one at each step, we cannot repeat this process indefinitely. Therefore, for some 
k, neither Case 1 nor Case 2 holds for S, and thus for all < j < k, Ej D -Ej+i and Fj C -^j+i- 
This implies that 5" is a path space. □ 



Proof of Theorem 3.4 Let Ai be the set of path spaces described in the theorem. We first show. 



by contradiction, that all path spaces in Ad are maximal. Suppose not. Then there exists some 
M = U^^qO, G M that is strictly contained in the path space 5' = U^^qO-. 

Suppose Oj C O'l for some < j < k and some < I' < k'. Since Ti and T2 have no 
splits in common, Ej D F^ = Fj f] E[ = 0. This implies that Ej C E'^ and Fj C F^. Also, 
F[ C Ct2{E[) ^ CtjC-Ej) = Fj, where the last equality follows from Condition 2 on the path spaces 
in Ai. Hence, F'^ = Fj. This implies E[ C Ct^{Fi) = CTi{Fj) = Ej, where the last equality follows 
from Condition 1. Therefore, E[ = Ej, and hence Oj = 0[. 

Therefore, no orthant in 5' can strictly contain an orthant from M, so S' is exactly the orthants 
forming M as well as at least one other orthant. Let j be the smallest index such that the orthant 
Oj-i is in M and 5', but O'j, O'j^^, Oj^i_^ are not in M and Oj = 0'j_^i. Since Oj is an orthant 
distinct from those in M, Ej^i = E'-^-^ D E'- D E'-^i = Ej and Fj^i = F'-_^ C Fj C F'._^^ = Fj. The 
crossing set in of the splits added as we transition from Oj-i to O'^ is contained in the set of 
splits dropped at this transition. That is, XEj_i{F'-\Fj-i) C Ej-i\Ej C Ej-i\Ej. But Conditions 
1 and 3 imply that the crossing set of any element / G Fj\Fj-i is exactly the splits dropped at 
the j-th step, or XEj_^{f) = Ej^i\Ej. In particular, for every / G Fj\Fj„i C Fj\Fj^i, we have 
Xej^xU) — Ej-i\Ej- This implies that XEj_-i{Fj\Fj^i) = Ej^i\Ej, which contradicts the strict 
inclusion that we just showed. So no path space in Ai is contained in another path space. 

Let S = U^^qCj be some path space that is not in Ai. We will now prove that 5 is contained 
in another path space, S', and hence is not maximal. Since 5 ^ A^, at least one of the three 
conditions does not hold. 

Case 1: There exists a < j < A; such that E' = CTi{Fj)\Ej is not empty. That is. Condition 
1 does not hold. 

We now construct a path space in which the splits E' are dropped at the j-th step instead of an 
earlier one. Define S' = uL.qO'^, where 



O' 



Oi + 0{E') if < i < j 
Oi '\ij<i<k 



Since Oi C 0[ for all i ^ j and Oj C Oj + 0{E'), we have S C S'. It remains to show that S' is 
a path space. By definition, E' is compatible with Fj D F'j-i ^ ••• ^ -^b, so E'- and F- = Fi are 
compatible for all < i < j. The orthants remain unchanged for j < i < k. Since Et^ = E'q D 
E'l 5 ... 5 Ej D ... D E'f^ = 2), S' is a relaxed path space, and hence a path space by Lemma 
Therefore, S is strictly contained in the path space S' . 

Case 2: There exists < j < k such that F' = CT2{Ej)\Fj is not empty. That is. Condition 2 
does not hold. 



3.6 
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We will now construct a path space in which the splits F' are added to the tree at the j-th step, 

O' 



instead of a later step. Define S' = U^^qO', where 



O' 




if < i < J 
if J < i < A: 



Since Oi C O'- for all i ^ j and Oj C Oj + 0{F'), we have S C S". It remains to show that 
S' is a path space. By definition, F' is compatible with Ej D ^j+i ^ ••• ^ ^fci so F^' and 
= E'j are compatible for all i > j. The orthants remained unchanged for < i < j. Since 
Fq C F{ C ... C Fj C F'-j^^ C ... C F^, 5' is a relaxed path space, and hence a path space by 



Lemma 3^ Thus, S is strictly contained in the path space 5". 

Case 3: Neither Case 1 nor Case 2 holds, and, for some 1 < j < A;, there exists / G Fj\Fj^i and 
a minimal element g in P{T{Ej^i),T{ET2\Fj-i)) such that g < f in P{T{Ej^i),T{ET2\Fj-i)). 
That is. Conditions 1 and 2 hold, but Condition 3 does not hold. We will now construct a path 
space in which we add the splits g and / in two distinct steps, instead of during the same step. 
Define S' = uf+(}0^, where 

'Oi if < i < j 

0'i=lo {{Ei.i\XE^,,{g)) U F_i U g) if i =j 

Oi-i ifj<i<A;-l 

We will first show that O'j is neither contained in nor contains any orthant from S. We must have 
^Ej-iig) 7^ 01 or else g S CTaC-Fj^ implying Case 2 holds, which is a contradiction. This 
implies that Ej^i D Ej^i\XE^_^{g), or E'-_^ D E'-. Since 5 < / in P{T{Ej^i),T{ET2\Fj-i)), we 
have that XEj^^{g) C XEj_iif)- To add / at step j, we must drop any splits in Fj_i that are 
incompatible with /, which implies XEj_^{f) ^ Ej-i\Ej. This, along with the previous statement, 
implies that C Ej^i\Ej. In turn, this implies that Ej C Ej-i\XE^_^{g), or Ej_^_i C Ej. 

Therefore, we have shown that Ej_^ D Ej D ^j+i^ &s desired. 

Since g £ Fj-aVFj-i, we have that Fj_^ = Fj-i C Fj_i U 51 = Fj, and hence it remains to 
show that Fj C Fj_^^. First, we will show that g G Fj. Since g < f, XEj_i{g) C XEj_i{f) C 
^E -x{Fj) ^ This implies that (7 S CT2{Ej) = Fj by Condition 2. Next we will show 

that Fj = F,_iU5 C F,- = Fj+i. For any /' G F,_i U g, C Xt,(Fj_i) U Xt^Is) C 

XxiiFj) by definition of closure and g G Fj. Along with the partitioning property, this implies that 
XTiif) n CtAFj) = 0, or Xt AD n Fj = by Condition 1. Then /' G Cr^iEj) = Fj, as desired. 
Furt hermore, Xe^_,(F~^) = Xe^^AFj-i) U ^E^^^g) = U Xe^_^) C and thus 

/ ^ Fj_i U g by definition of the closure operator. So / E Fj\Fj-i U g, and hence Fj C Fj_|_^. 
Therefore, Fj_i C Fj C F,+i. 

It remains to show that the splits in Oj are mutually compatible. By the definitions, 
CTi(Fj-iUff) = Cti(F,_i) n CTi(ff) 3 Fj_i\Xti(5) 2 Fj_i\Xi=;^._^(ff), and hence the splits of 
Oj are mutually compatible. The other orthants remain unchanged, and thus 5 is a path space. 
Therefore, S' is a path space that strictly contains S, so S is not maximal. □ 

Recall that in a poset F, x < y is a cover relation, or y covers x, if there does not exist any 
z €z P such that x < z < y. A chain is a totally ordered subset of a poset. A chain is maximal 
when no other elements from P can be added to that subset. See j251 Chapter 3] for an exposition 
of partially ordered sets. 

Theorem 3.7. Let g : K{Ti,T2) Tn he given by g{A) = Oa, where Oa = 0{CtAA) U A), for 
any A G K{Ti,T2). For any maximal chain Aq < Ai < ... < A^ in K{Ti,T2), define h{AQ < Ai < 
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... < Ak) = uf^o9(^i)- Then uf^gS'l^j) = U^^qO^. is a maximal path space and h is a bijection 
between maximal path spaces from Ti to T2 and maximal chains in K(Ti, T2). 

Proof. The map g is one-to-one, because A ^ A', then Xx^iA) ^ Xti{A'), and hence Cti{A) ^ 
Cti{A') by the partitioning property. We now show that h maps maximal chains in K{Ti,T2) to 
maximal path spaces. 

Let = Aq < Ai < ... < Ak = be a maximal chain in i^(Ti,T2). For every < i < /c, 
let Fi = Ai and Ei = CTi(^i). We now show that U^^qCj is a path space. Since K{Ti,T2) 
is the closed sets of £^7^2 ordered by inclusion, Ai C Ai^i for all < i < k. This implies that 
= Fq C Fi C ... C Fk = E{T2), since Fi = Ai for ah i. It also implies that Xri(^j) C 
If XxiiAi) = X-jPj (^j+i), then Ai^i C = Aj, by definition of the closure and since Ai is a closed 
set. This is a contradiction, and therefore, Xx^iAi) C ^^'^(jdj+i). This implies Ei = CxiiAi) D 
C7Ti(Aj+i) = ^Jj+i, and hence Eq D Ei D ... D -Efc. Also, Eq = Ct^^Aq) = Cti(0) = Et^, and 
-E'fc = CTi{Ak) = Cti{Et2) = since otherwise, T2 could contain more than n — 2 splits. Finally, 
for all < i < A;, Ei is compatible with Fi by definition. Therefore, ^'j- nOjE i U Fj) is a path space. 



We will now show that ujLgO,; satisfies the 3 conditions of Theorem 3.4 and hence is maximal. 



Since Ei = CT^{Fi), Condition 1 is met. As in any path space, Fi C CTziEi). For any / S CTiiEi), 
Xr^if) n Ei = 0. Since Xr^^^Ai) and Cr^iAi) partition Et^ and Ei = Cr^^Ai), then Xr^lf) C 
XxiiAi). So by definition f £ Ai = Ai = Fi. Therefore, Fi 5 CT2{Fi) also, and Condition 2 holds. 

To show Condition 3, suppose that for some 1 < i < A;, there exists / G Fj\Fj^i and a 
minimal element g in P{T(Ej^i), T{E'r2\Fj^i)) such that g < f in P{T{Ej^i), T{ET2\Fj-i)). Then 
XrAa) C Xti(/), so / ^ Fi-i U g. This implies = Fi_i < F~^ < Fi_i UgUf < Fi = Ai, 
and hence Ai < Ai^i is not a cover relation, which is a contradiction. Therefore, Condition 3 also 
holds, and ufLoOi is a maximal path space. 

So as claimed, if ylo < ^1 < ••• < ^fc is a maximal chain, then h{Ao < ... < A^) is a maximal 
path space. It remains to show that his a bijection. For any maximal path space U^^qCj, Fi < .Fj+i 
is a cover relation for all < i < A; since for any / G Fi^i\Fi, FiU f = -Fj+i by Condition 3 of 



Theorem 3.4 This implies that = i^o < -^1 < ••• < -^fc = E{T2) is a maximal chain in K(Ti,T2) 
such that h{FQ < Fi < ... < Ff^) = U^^qCj, and hence h is onto. We have that h is one-to-one, 
because g is one-to-one. Therefore, h is a bijection, which establishes the correspondence. □ 
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4" '6" 'n-2" 

'2-i<=l''=2! {62,63,64} {64,65,6^,} {6„_2, e„_3, 6„_4, 





fn-2-<V2^ 



(a) Tree Ti. (b) Tree T2. (c) Incompatibility poset P{Ti,T2). 

Figure 5: A family of trees whose path poset is exponential in the number of leaves. 

The number of elements in a path poset K{Ti, T2) can be exponential in the number of leaves 
in Ti and T2. For example, for any even positive integer n, consider the trees Ti and T2 depicted 



in Figures 5(a) and 5(b) Their incompatibility poset is given in Figure 5(c), Let W be the set of 



minimal elements in P(Ti,r2). Then \W\ = Each subset of is a distinct closed set, and 

hence an element in K(Ti,T2). This implies there are at least 2^~ elements in K{Ti,T2), and 
hence also an exponential number of maximal chains. 
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4 Geodesies in Path Spaces 



Given a path space, this section shows how to find the locahy shortest path, or path space geodesic, 
between Ti and T2 within that space in hnear time. We do this by transforming the problem into 



a Euclidean shortest-path problem with obstacles ([T7] and references) in Theorem 4.4 We next 



reformulate the problem as a touring problem [91. A touring problem asks for the shortest path 



through Euclidean space that visits a sequence of regions in the prescribed order. Lemmas 4.8 and 



4.9 give conditions on the path solving the touring problem. The linear algorithm for computing 



the path space geodesic is given in Section 4.2.1, with Theorem 4.10 proving its correctness. 



4.1 Two Equivalent Euclidean Space Problems 

Let Ti and T2 be two trees with no common splits, and let S = U^^QO{Ei U Fi) be a path space 
between them. Define the path space geodesic between Ti and T2 through S to be the shortest path 
between Ti and T2 contained in S. Let ds{Ti,T2) be the length of this path. 

We will now show that the path space geodesic between Ti and T2 through a path space 
containing A; + 1 orthants is contained in a subspace of 7^ isometric to a subset of a lower or equal 
dimension Euclidean space, V{R''). For <i<k, define the orthant 



Vi = {{xi, Xfc) G R'' : Xj < j < i and Xj > if j > i}. 



Let V{R^) = UjLo^i- 



We prove three properties of path space geodesies, and hence also geodesies, in Proposition 4.1 



Proposition 4.2 and Corollary |4.3| These properties will imply that if a set splits shrink to or start 
growing from length at the same point on the path space geodesic, then we know the length of 
each edge at any other point on the path space geodesic. Analogous properties were proven by 
Vogtmann [2^ for geodesies. 

Proposition 4.1. The path space geodesic is a straight line in each orthant that it traverses. 

Proof. If not, replace the path within each orthant with a straight line, which enters and exits the 
orthant at the same points as the original path, to get a shorter path. □ 

Proposition 4.2. Moving along the path space geodesic, the length of each non-zero edge changes 
in the trees on it at a constant rate with respect to the geodesic arc length. 



Proof. By Proposition 4.1 each edge must shrink or grow at a constant rate with respect to the 
other edges within each orthant, but these rates can differ between orthants. So it suffices to 
consider when the geodesic goes through the interiors of the two adjacent orthants Oi = 0{EiU Fi) 
and Oj+i = ©(E'j+i U-Fj+i), and bends in the intersection of these two orthants. Let a be the point 
at which the geodesic enters Oi, and let b be the point at which the geodesic leaves Oi+i. 

The edges Si\£'j+i are dropped and the edges Fi+i\Fi are added as the geodesic moves from 
Oi to Oj+i. Thus the edges £'i\-Ej+i and all have length in the intersection 0{Ei+iU Fi). 

Let m = I Ei+i UFi\, the dimension of Oi n Oi+i . Consider the subset S = Ha U Hb of Oi U Oj+i , 
where Ha is the affine hull of a U {Oi Ci Oj+i) intersected with Oi and Hb is the affine hull of 
b U {Oi n Oi+i) intersected with Oj+i. This subset can be isometrically mapped into two orthants 
in M'"^-'^ as follows. For each tree T G Ha, let the first m coordinates be given by the projection of 
T onto Oi n Oi+i. Let the (m + l)-st coordinate be the length of the projection of T orthogonal to 
OiOOi+i. More specifically, let the edg es in Ei^i U Fi be ei , 62, • . em- Then we map T to the point 
(IciIt, |e2|T) ••••) |em|r, s) in M™"*"^, where s = \/Yle&Ei\Ei+i\^\T- Similarly, for each tree T G Hb, 
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let the first m coordinates be given by the projection of T onto Oi n Oj+i. Let the (m + l)-st 
coordinate be the negative of the length of the projection of T orthogonal to Oi n Oi+i- In other 



words, we map T to the point (|ei|T, \^2\t-, \em\T, —s) in M"^+ , where s = y X]eeFi+i\Fi I'^iT- 

We have mapped 5 into Euclidean space, and hence the shortest path between the image of a 
and the image of b is the straight line between them. Along this line, each edge ei, changes 
at the same rate with respect to the geodesic arc length. Since we can make this argument for each 
pair of consecutive orthants, we have proven this proposition. □ 

Corollary 4.3. Let T be a tree on the path space geodesic between Ti and T2 through the path 



space Uf^QO{Ei U Fi 
/i,/2 G Fj\Fj^i, and if i < j < k, we have 



Suppose T £ Oi. Then if 1 < j < i, we have jj^j^ 



\ei\T 



\e2\T 

\e2 



1/2 h 
1/2 It2 



for any 



jl^ for any 61,62 G Ej^i\Ej. 



Proof. Proposition 4.2 implies that the length of each edge in Ti shrinks at a constant rate until it 
reaches as we travel along the path space geodesic, and the length of each edge in T2 grows at a 
constant rate from starting at some point along the path space geodesic. Since for any 1 < j < k, 
the edges Ej-i\Ej reach length at the same point along the path space geodesic, each edge in 
Ej^i\Ej must be changing at a constant rate with respect to the lengths of the other edges in 
Ej^i\Ej. Similarly, since the edges Fj\Fj^i start growing from at the same point along the path 
space geodesic, each edge in Fj\Fj-i is changing at a constant rate with respect to the lengths of 
the other edges in Fj\Fj^i. □ 

Therefore, there is one degree of freedom for each set of edges dropped, or alternatively for 
each set of edges added, at the transition between orthants. Thus, the path space geodesic lies 
in a space of dimension equal to the number of transitions between orthants. We will now show 
that each path space geodesic lives in a space isometric to V(R'^). For example, in Figure 6(a) , the 
path space Q consists of the orthants 0(161, 62, 63}), ^^({/i) £2, 63}), and 0{{fi, f2, fs})- We apply 
Theorem |4.4| to see that the geodesic through Q is contained in the shaded region of shown in 
Figure |6(b)[ 



(61,62,63) 





geodesic 

(a) Part of T5. (b) Isometric mapping to V(IR^) 

Figure 6: An isometric map between a path space and y(M^). 



Theorem 4.4. Let Q = U*LoC'(^i U Fj) be a path space between Ti and T2, two trees in Tn with no 
common splits. Then the path space geodesic between Ti and T2 through Q is contained in a space 
isometric to y(]R*^). 
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Proof. For all 1 < j < k, let Aj = Ej_i\Ej and let Bj = Fj\Fj^i. Then {^jjLi and {Bj} 



are partitions of Et^ and Et^ respectively. By Corollary 4.3, any tree T' G Q on the path space 
geodesic satisfies the following property for each 1 < j < k: 



1. if T' E Oj and j < i, then there exists a Cj = Cj(T') > 0, depending on T' , such that 
for ah f e Bj, 



It' 



Ci 



I/It, - "J 



2. if T' G Oi and j > i, then there exists a = dj{T') > 0, depending on T', such that j^j^ = dj 
for all e G Aj. 

Let Q' C 7^ be the set of trees satisfying this property. For < i < n, define hi : Q' f] Oi ^ Vi hy 

h^{T') = hi{T{ci ■ Bi U ... Uc^-BiU di+i ■ Ai+i U ... U 4 • At)) 
= (-ci||Bi||, -Ci\\Bi\\, di+i\\Ai+i\\, dk\\Ak\\) . 

We claim that hi is a bijection from Q' n Oi to the orthant Vi in y(M^). The orthant Oj contains 
trees with = \Bi\ + \B2\ + ... + \Bi\ + |Ai+i| + .... + \Ak\ edges, and hence is an A^-dimensional 
orthant. All trees in Oi contain exactly the edges {Bi, Bi, Ai+i, A^}, so without loss of 
generality we can assign each edge to a coordinate axes so that the edges in Bi are assigned to 
coordinates 1 to | -Bi | , the edges in B2 are assigned to coordinates | -Bi | + 1 to | i?i | + | i?2 1 1 the edges 
in Ai^i are assigned to the coordinates + |i?2| + ••• + |-Bj| + 1 to + |i32| + .•• + + |Ai+i|, 
etc. Let Bj be the edge assigned to the j-th coordinate. By abuse of notation, for all 1 < j < i, 
let Bj be the A^-dimensional vector with a in every coordinate except those corresponding to the 
edges in Bj, which contain the lengths of those edges in T2. Similarly, for all i < j < k, let Aj be 
the A^-dimensional vector with a in every coordinate except those corresponding to the edges in 
Aj, which contain the lengths of those edges in Ti. For example, Bi is the A^-dimensional vector 
(l/i|r2, 1/2IT2, I/|Bi||t2,0, ...,0). 

Then Q' n Oi is generated by the vectors | , UBfu ' ' ' " ' ' ||a'+i|| ' • • ' ||Ak|| | • Since these gen- 
erating vectors are pairwise orthogonal, they are independent, and hence Q' nOi is a A;-dimensional 
orthant contained in Oi. Furthermore, for all 1 < j < i, corresponds to the tree T (^yp7]| " Bj 

and for all i < j < k, j^^j^ corresponds to the tree T ' ^i) " ^ ^ j ^ k, let Uj be the 

fe-dimensional unit vector with a 1 in the j-th coordinate. Then for 1 < j < i. 

Similarly, for all i < j < k. 

The basis of Vi is {— ui, — Uj, Ui+i, Uk}, so hi maps each basis element of Q' n to a unique 
basis element of Vi. Thus, hi is a linear transformation, whose corresponding matrix is the identity 
matrix, and hence a bijection between Q' n Qi and Vi for all i. Furthermore, since the determinant 
of the matrix of hi is 1, hi is also an isometry. So Q' is piecewise linearly isometric to V{M.^). 

For all < i < n, the inverse of hi is gi : Vi ^ Q' defined by gi{—xi,... — Xi, Xi+i, Xk) = T' , 
where Xj > for all 1 < j < A; and T is the tree with edges Ei U Fi with lengths • \e\T2 if 

\x\ I I . ^ 

e E Bj for 1 < j < i and TT^TiT • \e\Ti if e E Aj for i < j < k. 
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Notice that if T' e Q' D Oi D Oi+i, then hi{T') = hi+i{T'), since the lengths of all the edges 
in Ai+i and Bi+i are 0. Therefore, define h : Q' ^ V{R^) to be h{T') = hi{T') if T' e n Q', 
which is well-defined. Define g : V{'R'') — > Q' by setting g{—xi, ... — Xi, Xj+i, Xk) = gi{—xi, ... — 
Xi, Xi+i, Xfc), for all 1 < i < k and for all Xj > for all 1 < j < A;. Then g is also well-defined 
and the inverse of h. 

For any geodesic q in Q', map it into VCR.^) by applying h to each point on q to get path p. 
Notice that since both hi and gi are distance preserving, p is the same length as q. We claim p 
is a geodesic in V(R^). To prove this, suppose not. Let p' be the geodesic in V{M.^) between the 
same endpoints as path p. Then p' is strictly shorter than p. Use g to map p' back to Q' to get 
q' . Again distance is preserved, so q' is strictly shorter than q. But q was a geodesic, and hence 
the shortest path between those two endpoints in Q' , so we have a contradiction. Therefore, the 
geodesic between Ti and T2 in Q is isometric to the geodesic between A = (||Ai||, ||Ak||) and 
i3 = (-||Bi||,...,-||Bk||) in V{R^). □ 

Thus, finding the shortest path between two trees through a (A; -|- 1)— orthant path space is 
equivalent to finding the shortest path between a point A in the positive orthant and a point 
B in the negative orthant of 1/(M'^). This problem can be transformed into an obstacle-avoiding 
Euclidean shortest path problem by letting A and B be points in M*^, and letting the orthants which 
are not in V{M.'') be obstacles. We now formulate this problem as a touring problem. Let Pi be 
the boundary between the i-th. and {i + l)-st orthants in V{M.^), for all 1 < i < A;. That is. 

Pi = {{xi, Xk) G M'^ : Xj < if j < i; xj = if j = i; xj > if j > i}. 

Then our problem is equivalent to finding the shortest path between A and B in M'^ that intersects 
Pi, P2, Pk in that order. 

In dimensions 3 and higher, the Euclidean shortest path problem with obstacles is NP-hard in 
general [6], including when the obstacles are disjoint axis-aligned boxes [18]. The touring problem 
can be solved in polynomial time as a second order cone problem when the regions are polyhedra 
^20j. In the special case of the above touring problem, we find a simple linear algorithm. 



4.2 Touring Problem Solution 



First, Lemma 4.5 establishes when AB is the solution to our touring problem. Next we introduce 
the idea of a locally shortest, ordered path, and prove two conditions that all such paths must 



satisfy in Lemmas 4.8 and 4.9 Theorem 4.10 shows that repeatedly applying the second condition 
gives the linear algorithm for finding the shortest, ordered path from A to B. Throughout this 
section A = (ai, 02, a/c) and B = (— 61, — 62, — ^fc) will be points in M*^ with ai,bi > for all 
1 < i < A;. 

Lemma 4.5. The line AB passes through the regions Pi,P2, -..^Pk in that order and has distance 
'^=\lY.'l=M + h? if and only ^/ |i < g < ... < ^. 

Proof. Parametrize the line AB with respect to the variable i, so that t = at j4 and t = 1 at B, 
to get (xi, Xfc) = (ai, ...,ak) + t{—ai — 61, —ak — bk). Let ti be the value of t at the intersection 
of AB and Pi. Setting Xi = 0, and solving for t gives ti = „ °?h ■ For AB to cross Pi, P2, .... Pk in 
that order, we need ti < to < ... < tk or < < ... < . Since for any 1 < i, j < k, 

< a i-h- ^® equivalent to p < by cross multiplication, we get the desired condition. By the 

Euclidean distance formula, AB = A/X^iLi(ai -|- □ 
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Corollary 4.6. Let A = (ai, o^) and B = (—61, —hk) he points in M with ai, bi > for all 



1 < i < k. Then ^ = if and only if AB intersects Pi Pi Pj+i. 



Proof. This follows directly from the proof of Lemma 4.5 □ 



In general, we will not have |^ < < ... < and hence the shortest path is not a straight line. 
We now show how to find where a path bends, using the idea of locally shortest, ordered paths. 
These bends can be straightened by isometrically mapping the problem to a lower dimensional 
space, until Lemma |4.5| applies. 

An e-neighbourhood of a path is all points in M*^ within e > of at least one point on that path. 
A locally shortest, ordered path, q, is a path from A to B which passes through Pi,...,Pk in that 
order, and for which there exists some e > such that there is no shorter path q' from A to B 
contained in the e-neighbourhood of q that also passes through Pi, in that order. For all i, let 
Pi be the first point at which the locally shortest, ordered path under consideration intersects Pi. 
Then it is easy to show that for all 1 < i < k — 1, any locally shortest, ordered path g is a straight 
line, possibly of length 0, between pi and Pi+i, and q intersects each Pi at exactly one point, pi. 



The following corollary of Theorem 4.4 explains when a touring problem can be isometrically 
mapped to a lower dimension. 

Corollary 4.7. Consider a locally shortest, ordered path from A = (01,02, ..-,0*;) to B = (—61, —62; 
...,—bk) through Pi, Pk. Let {Mj}^-^ be any ordered partition 0/ {1, 2, fe} such that i, I G Mj 
implies pi = pi. Then this path is contained in a region ofM.'^ isometric to V{W^). 

Proof. Suppose i,i + 1 are in the same block in {Mj}JL^. Then pi = pj+i, and travelling along 
the pre-image of the path in tree space, the tree loses splits Ei-i\Ei and Ei\Eij^i simultaneously, 
and gains splits Fi\Fi^i and Fj+i\Pj simultaneousl y. H ence, this path is in the path space Oq U 

Uf=iO {{nieM^E,) U {UieM.Fi))). Apply Theorem 



4.4 



to get the desired result. □ 



Notice that under the mapping to V(R"^) described in the above proof, A is 



mapped to .4 = (^EieA/i , -^EieMa , -^EieM™ a?j and B is mapped to B = 

The following two lemmas give simple constraints on locally shortest, ordered paths. 

Lemma 4.8. Let q be a locally shortest path from A to B that passes through Pi, P2, Pk in that 
order. Let pj be the intersection of q and Pj for each 1 < j < k. If j^j < 1^ ■■■ < for some 
1 < J < i < k, q is a straight line until it bends at pj = pj+i = ... = pi, and pj-i 7^ pj, then 
Pi = Pi+i- 

Proof. This proof is by contradiction, so assume that pi 7^ Pi+i- By properties of locally shortest, 
ordered paths, g is a straight line from pi to Pi+i- Let Y = (— yi, —yi, Vi+i, ■■■,yk)i where yj > 
for all 1 < j < /c, be a point on the line pipi+i, e > past pi. We will now show that AY intersects 
Pi, P2, Pi in that order. 

Parametrize the paths q and AY with respect to time t, so that t = at A and t = 1 at Y. 
The j-th coordinate, for 1 < j < J — 1, decreases linearly from aj to —yj in both q and AY, and 
thus become at the same time in both paths. This implies that since q crosses Pi, Pj_i in that 
order, AY also crosses Pi, ...,Pj_i in that order. 

Let tj be the time at which AY intersects Pj, for 1 < j < i. Then = Oj + tj{—yj — aj) or 
ti = — — • In each coordinate between J and i becomes at the same time. These coordinates 
then decrease linearly, so the ratio between any two consecutive coordinates remains constant as 
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time increases. This implies = for each J < j < i. Since < < ... < ^ hy 

Vj+i "j+i ~ bj — bj+i — — bi 

the hypothesis, then < < ... < This implies < — — < ... < or 



y.7 — y.J +1 — "' — Vi' ^ a.j+yj — ajj^i+yjj^i — — ai+y, 

tj < < ... < ii- Thus AY intersects Pj, Pj+i, Pi in that order. 

It remains to show that AY intersects Pj-i before Pj, which we do by contradiction. So assume 
that tj < tj-i. Let rj_i and rj be the points of intersection of AY with Pj-i and Pj, respectively. 
By the hypotheses and assumption, rj and pj are contained in Pj\Pj^i. Since Pj~i and Pj are 
convex, rJZipJ^ and rjpj are contained in Pj-i and Pj, respectively. Now rJZipJZi intersects 
rjpj inside the triangle AYpj. This implies that rJpJ passes from Pj\Pj-i into Pj-i DPj, on the 
boundary of Pj, and back into Pj\Pj-i. But this contradicts the convexity of Pj. Thus < tj, 
and AY passes through Pi, P2, Pi in that order. 

By the triangle inequality, AY is shorter than the section of q from A to Y. This contradicts q 
being a locally shortest, ordered path, and thus pi = Pi+i. □ 

Lemma 4.9. For any locally shortest paths from A to B that pass through Pi, P2, Pk in that 
order, «/|^<|^<...<^> then any such path intersects Pi Ci Pi+i. 

Proof. Suppose that ^ < ^ < ... < ^ and consider some locally shortest path, q, such that 
Pi 7^ Pi+i. Parametrize q with respect to the variable t, so that the path starts at A when t = 0, 
ends at B when t = 1, and passes through Pj at point pj = {pj,i,Pj,2, ■■■,Pj,k) when t = tj, for 
all I < j < k. If Pj = Pj+i for some 1 < j < i and q bends at this point, then by repeated 



applications of Lemma 4.8 it also passes through Pi n P^+i and we are done. So assume that q is 
a straight line from A to Pi+i. Thus, the i-th coordinate changes linearly from aj to —bi, and from 
the parametrization of this, we get tj+i = ^t^'' • 

Case 1: 7^ (That is, the locally shortest ordered path does not bend at Pi+i.) 

In this case, = = Oj+i + - Oj+i), which implies tj+i = ^i+i+bi+i ■ Equate 

this value of tj+i with the one found above, and rearrange to get Pi+i^i = Oj — The 
definition of Pj+i and the assumption pi ^ pj+i implies that pi+n < 0. Hence, < !h±li^h+hl^ 
which can be rearranged to < a contradiction. 

Case 2: pj+i_i+2 = (That is, the locally shortest ordered path bends at Pi+i, and pj+i = Pi+2-) 

Let J > 2 be the largest integer such that pi+j = Pi+i, but pj+j+i 7^ Pi+i. Apply Corollary 4.7 
using the partition {1}, {2}, {i}, {i + 1}, {i + 2, + J},{i + J + 1}, {k} to^reducejhe space 
by J — 2 dimensions. A and B are mapped to ^ = (oi, afc_.(j_2)) and B = {—bi, — 6fc_(j_2)), 
respectively, in the lower dimension space, where: 



if j < i + 1 
12i=2(4+i ifj = i + 2 and bj 
aj+.j-2 ifj>i + 2 



bj if j < i + 1 

T.I2 bli ifj = i + 2 



+1 

bj+j-2 if j > i + 2 



Let k = k—{J—2). Let pj be the image oipj in M'^ under the above mapping if j < i+2 and the image 
of Pj+j-2 li 3 > i + 2. Let Pj = {(xi, ■■■,x~} e R'^ : xi < ifl < j; xi = if I = j; xi > if I > j}. 

So Pj is the boundary between the j-th. and {j + l)-st orthants in the lower dimension space M.^. 
Let q be the image of q. 



Then g is a straight line from A to Pi+i, and pj+i = pi+2 / Pi+3, so q 



Since q does not intersect Pj_|_2 H -Pj+3, by the contrapositive of Lemma 



4.8 



sends in Pj+i Pi Pi+2- 
' > In R^, 
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this translates into the condition that > J ' =. Cross-multiply, square each side, add 

a^+ifti+i, and rearrange to get ^ > \ • 

The remaining analysis is in M^. If the locally shortest, ordered path is a straight line through 
Pi+i, then we make the same argument as in Case 1. Otherwise, since the path does not bend 
at Pi, the i-th coordinate changes linearly from to —hi. We use this parametrization to find 

ti+2 - ti+i - ^,7+5—- 

Furthermore, the {i + l)-st to {i + J)-th coordinates decrease at the same rate from A to Pi+i 
and at the same, but possibly different than the first, rate from pi+i to B. Therefore, we can 



apply Corollary 4.7 to the partition {1}, {2}, {i}, {i + 1, i + 2, z + J}, {i + J + 1}, {/c} to 



isometrically map the locally shortest, ordered path into R™ ^\ Let a = y'^i=iCii^i, and let 

b = \l ■ Then in W^~^'^~^\ the {i + l)-st coordinate of the locally shortest ordered 
path changes at a constant rate from a to —b. This implies = a + 6 — a), or tj+i = 

Equate the two expressions for tj+i to get pi+i i = ai — ^"'J" • By definition of -Pj+i, Pi+i i < 0. 

This implies ^ < £ = y^'-^"'+' . But we showed that y^'-^"'+' < ^ go ^ < 
which is also a contradiction. □ 

By repeatedly applying this lemma, we find the lower dimensional space that all the locally 
shortest, ordered paths lie in. In this space, the ratios derived from the coordinates of the images 
of A and B form a non-descending sequence. The following theorem gives the shortest path through 
ViM}') from a point in the positive orthant to a point in the negative orthant, or equivalently, the 
shortest tour that passes through Pi, in M'^. 

Theorem 4.10. Let A = (ai, 02, Ofc) and B = (—61, —h2, —bk) with Ui, 5j > for alll < i < k 
be points in M*^. Alternate between applying Lemma 4-9 and Corollary 4.I until there is a non- 



descending sequence of ratios < ^ < ... < if^ , where at and bi are the coordinates in the lower 

bl b2 bm ^ ^ ^ ^ 

dimensional space. There is a unique shortest path between A = (oi, ....,am) and B = (—61, —bm) 



in y(M'"), with distance yX^ilLiC'^i + bi)^. This is the length of the shortest ordered path between 
A and B in V{R''). 



Proof. For the smallest i such that > Lemma 4.9 implies that pi = pi+i in all locally 

shortest, ordered paths in M'^. Thus, we can apply Corollary 4.7 using the partition {l},{2},...,{i — 
1}, {i, i + l}, {i + 2}, {m} to reduce the space containing all locally shortest, ordered paths by one 
dimension. Repeat the previous two steps until the ratio sequence in the lower dime nsional space 



4.5 



the geodesic 



is non-descending. Let < < ... < 3?^ he this ratio sequence. By Lemma 

61 - 62 ~ ~ 

between A and i? is a straight line, and hence unique. Furthermore, its length is y'SI^iC'^* + ^«)^- 
Since we mapped from V{R^) to V{R"^) by repeated isometries, both the length of the path and 
the order it passes through Pi,...,Pm, or their images, remain the same. The straight line is the 
only locally shortest, ordered path in M™', so its pre-image is the only locally shortest, ordered path 
in M'^ and thus must be the globally shortest path. □ 
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4.2.1 PathSpaceGeo: A Linear Algorithm for Computing Path Space Geodesies 



Theorem 4.10| can be translated into a linear algorithm called PathSpaceGeo, for computing the 



path space geodesic between Ti and T2 through some path space S = uf-^oOfc. For all 1 < i < /c, 
let Ai = Ei^i\Ei and let Bi = Fi\Fi-i. Let ai = \\Ai\\ and bi = \\Bi\\, for all 1 < z < A;. 
Let 1 < i < A; be the least integer such that ^ > Then by Theorem 



1 



4.10 



to find 



7 To the ra- 



the path space geodesic through S, we should apply Lemma 4.9 and Corollary 
tio sequence to map the problem to ^(M'^^"'^), where the ratio sequence becomes 

1^, ^j^2_)_fc2^^ "' F^f T^- Repeat this process until the ratio sequence is non-descending. 

Unfortunately, this process is not deterministic, in that different non-descending ratio sequences 
can be found for the same geodesic, depending on the starting path space. This occurs, because 



by Corollary 4.6, two equal ratios can be combined to give a ratio sequence corresponding to a 
path with the same length. However, if we modify the algorithm to also combine equal ratios, the 
output ascending ratio sequence will be unique for a given geodesic. 

Alternatively, define the carrier of the path space geodesic through S between Ti and T2 to be 
the path space Q = U^^QOc(i) ^ ^ such that the path space geodesic through S traverses the 
relative interiors of 0^(0), C'c(i)) •••) ^c{i)i where the function c : {0, 1, /} — > {0, k} takes i to 
c(i) if the i-th orthant is Q is the c(i)-th orthant in S. If a path space geodesic is the geodesic, we 
just write carrier of the geodesic. The carrier of the path space geodesic is the path space whose 
corresponding ratio sequence is the unique ascending ratio sequence for the path space geodesic. 

We now explicitly describe the algorithm for computing the ascending ratio sequence corre- 
sponding to the path space geodesic, PathSpaceGeo, and prove it has linear runtime. 

PathSpaceGeo 

Input: Path space S or its corresponding ratio sequence R = ^, 

Output: The path space geodesic, represented as an ascending ratio sequence, which is understood 
to be the partition of R where the ratio V^£=LJli corresponds to the block l^, Tr^\. 

Algorithm: Starting with the ratio pair l^^^, ff }' PathSpaceGeo compares consecutive ratios. If 

for the i-th pair, we have > then combine the two ratios by replacing them by /^=^= 

in the ratio sequence. Compare this new, combined ratio with the previous ratio in the sequence, 
and combine these two ratios if they are not ascending. Again the newly combined ratio must be 
compared with the ratio before it in the sequence, and so on. Once the last combined ratio is 
strictly greater then the previous one in the sequence, we again start moving forward through the 
ratio sequence, comparing consecutive ratios. The algorithm ends when it reaches the end of the 
ratio sequence, and the ratios form an ascending ratio sequence. 

Theorem 4.11. PathSpaceGeo has complexity @{k), where k + 1 is the number of orthants in 
the path space between Ti and T2. 

Proof. We first show the complexity is 0{k). Combining two ratios reduces the number of ratios 
by 1, so this operation is done at most k — 1 = 0{k) times. It remains to count the number of 
comparisons between ratios. Each ratio is involved in a comparison when it is first encountered 
in the sequence. There are k — 1 such comparisons. All other comparisons occur after ratios 
are combined, so there are at most /c — 1 of these comparisons. Therefore, PathSpaceGeo has 
complexity 0{k). Any algorithm must make k — \ comparisons to ensure the ratios are in ascending 
order, so the complexity is r2(A;), and thus this bound is tight. □ 
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We use PathSpaceGeo(S) = PathSpaceGeo^^, ^, to represent the ascending ratio 
sequence output by running PathSpaceGeo with input S. When PathSpaceGeo is run with an 
ascending ratio sequence as input, the output is the same sequence. Thus 



PathSpaceGeo ^PathSpaceGeo 
= PathSpaceGeo 



.....^...^^^^^ V^i' "' bj J ' ^j+i' "' bk 

Ol 02 Ok 

PathSpaceGeo lets us quickly calculate the shortest path through a maximal path space. 



5 Algorithms 



In this section, we show in Theorem 5.1 how to compute the geodesic distance between two trees 
Ti and T2 by computing the geodesic between certain smaller, related trees. This allows us to use 
the results from Sections 3 and 4, as well as either dynamic programming or divide and conquer 
techniques, to devise two algorithms for finding the geodesic between two trees with no common 
splits. Experiments on random trees show these algorithms are exponential, but practical on trees 
with up to 40 leaves and significantly better than the only other algorithm [15 published, to my 
knowledge. We also apply these algorithms to some biological data. 

5.1 A Relation between Geodesies 

Let Ti and T2 be two trees in 7^, with no common splits. The following theorem shows that there 
exists a maximal path space M containing the geodesic between Ti and T2 such that a certain 
subspace of it contains the geodesic between two smaller, related trees, T[ and Tg. As T[ and T2 
have fewer splits than Ti and T2, it is easier to compute this geodesic. Therefore, we can find the 
geodesic between Ti and T2, by finding the geodesic between all such T[ and Tg. 

Theorem 5.1. Let Ti and T2 be two trees in Tn with no common splits. Then there exist a 
maximal path space M = U^^gOj = ujLgC'(-Ej U Fi) between Ti and T2 and a maximal path space 
M' = ufjo O- = u'l~^0{Ei\Ek^i U Fi) between T[ = r(^Ti\^fc-i) and = T{Fk-i), such that 
M contains the geodesic between Ti and T2 and M' contains the geodesic between T[ and Tg. 

To prove this theorem, we first prove two lemmas which hold for any maximal path space 
M = uf^oOi = u'y^QO{Ei U Fi) between Ti and T2. In this context, let T[ = T{ET,\Ek_i) and 
T2 = T(Ffc_i). That is, T[ and T2 are exactly the trees Ti and T2 with the edges Ek-i and Fk\Fk-i 
contracted. Let M' = U^tTo = U^lQO{Ei\Ek_i U Fi) be the path space containing exactly the 
trees in M, with the edges Ek-i or Fk\Fk-i contracted. 



Lemma 5.2 shows that the carrier of the path space geodesic through M is contained in the 



orthants corresponding to the carrier of the path space geodesic through M' and Ok- Lemma 5.3 
shows that if M' does not contain the geodesic between T{ and Tg, and hence we can find another 
path space P' containing a shorter path space geodesic, then the corresponding path space P 
between Ti and T2 contains a path space geodesic at least as short as that through M. 

Lemma 5.2. Let Ti, T2, T[, T2, M and M' be as described above. Let Q' = U-^qO^^.^ be the carrier 
of the path space geodesic through M' . Let Q = {Q' (BO{Ek^i))UOk. Then ^^(Ti, T2) = dM{Ti,T2). 
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Proof. We want to show that PathSpaceGeo((5) = PathSpaceGeo(M). 
By properties of PathSpaceGeo, 



PathSpaceGeo(Q) = PathSpaceGeo ( PathSpaceGeo(Q')> Ittvi^^^ ) • 

V \\Fk\Fk-i\\ J 

Similarly, PathSpaceGeo(M) = PathSpaceGeo (^PathSpaceGeo(M'), ||r^^^) . But 

PathSpaceGeo(M') = PathSpaceGeo((5') by definition of the carrier of the path space geodesic, 
and thus PathSpaceGeo((5) = PathSpaceGeo(M). □ 

Lemma 5.3. Let Ti, T2, T[, Tg, M and M' be as described above. If M' does not contain the 
geodesic between T[ and in Tn, then there exists a path space P' between T[ and T2 such that 
dp,{T[,T^) < dM'{T[,T^) and dp{Ti,T2) < dM{Ti,T2), where P = {P' (B O(Sfe-i)) U Ok. 

Proof. Let Q' = u\^qO'^^-^ be the carrier of the path space geodesic through M'. Let q be the path 
space geodesic through Q' between T[ and Tg, and let qi = 0'^^-_^^ n Cc(j) H q for every 1 < i < I. 
Since q is not the geodesic from T[ to T2, q cannot be locally shortest in 7^. By Proposition 



4.1 



for all 1 < i < ? — 1, the part of q between q^ and gj+i is a line, and cannot be made shorter in 7^. 
Thus we can only find a locally shorter path in 7^ by varying q is the neighbourhood of some qj . In 
particular, there exists some e such that if s and t are the points on q, e before and after qj in the 
orthants Oc(j-i) ^-iid Oc(j)) respectively, then the geodesic between s and t does not follow q. Replace 
the part of q between s and t with the true geodesic between s and t to get a shorter path in 7^, with 
distance 4. Let C/ = 0{E'{uF^), O'^ = 0{E'^UF^), Oc(j) be the sequence of orthants 

through whose relative interiors the geodesic between s and t passes. Note that O", O'^ are not 
in M' . These orthants must form a path space, and thus P' = Q' U (U^^gO") is a path space. Since 
the path space geodesic is the shortest path through a path space, dpiiT^^Tl^) < dg < dQi{T[,T2). 
By definition of Q', dQi{T[,T2) = (ijv-/'(r(, T2), and hence dpi{T[,T2) < (ijv/'(T-[, T2), as desired. 

We will now show that dp{Ti,T2) < dM{Ti,T2). Let Q = (Q' © 0{Ek-i)) U Ok- Then Q C P, 
which implies dp{Ti,T2) < dQ{Ti,T2). By Lemma [5l2j dQ{Ti,T2) = dM{Ti,T2), and so P' is the 
desired path space. □ 



We use Lemma 5.3 to prove Theorem 5.1 



Proof of Theorem 5.1. Let S be any maximal path space containing the geodesic between Ti and 
T2. Suppose S consists of Z + 1 orthants, and let E be the set of edges dropped at the transition to 
the {I + l)-st orthant Oi. Let S' be defined by {S' © 0{E)) UOi = S. Then S' is a maximal path 



space between T[ and T2, as the conditions in Theorem 3.4 still hold. 

If S' contains the geodesic between T[ and Tj, then we are done. If not, then by Lemma |5.3[ 
there exists another path space P from Ti to T2, with dp{Ti,T2) < ds{Ti,T2) and dp'{T[,T2) < 
ds'{T[,T2). If P' does not contain the geodesic between T[ and T2, then we can keep applying 
Lemma |5.3| to a maximal path space containing it, until it does. As this process produces a path 
space containing a strictly shorter path space geodesic at each iteration, and there are only a 
finite number of path spaces between T[ and T2, it eventually finds the path space, Q' , containing 



the geodesic from T[ to Furthermore, letting Q = {Q' © 0{E)) U Oi, Lemma 5.3 implies that 
(iQ(Ti, T2) < ds{Ti,T2). But 5 contains the geodesic between Ti and T2, so ds{Ti,T2) < c?q(Ti, T2), 
and hence dQ(Ti,T2) = ds{Ti,T2). Let M' = U^~qO^ be a maximal path space between T[ and T2 
containing Q'. Then M = {M' © 0{E)) U Oi contains Q and is a maximal path space between Ti 
and T2 with the required properties. □ 
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We will now present two specific algorithms for computing geodesies. Each of these algorithms 



uses Theorem 5.1 to avoid computing the path space geodesic for each maximal path space between 
Ti and T2. This significantly decreases the run time. We call these algorithms GeodeMaps, 
which stands for GEOdesic DistancE via MAximal Path Spaces. The first algorithm uses dynamic 
programming techniques, and is denoted GeodeMaps-Dynamic, while the second uses a divide 
and conquer strategy, and is denoted GeodeMaps-Divide. 

5.2 GeodeMaps-Dynamic: a Dynamic Programming Algorithm 



Theorem |5.1| implies that for any element A in K{Ti,T2), we can find the geodesic between 
T{Xtj^{A)) and T{A) by only considering the maximal path spaces that contain a geodesic be- 
tween T{Xxj^{B)) and T{B) for some B covered by A. Furthermore, since given a geodesic qb 
between T(Xtj^(B)) and T(B) for some B covered by A, the candidate for the geodesic between 
T(Xtj^(A)) and T(A) is computed using only the carrier of gs, which is independent of the choice 
of maximal path space containing gs- 

This suggests that we can compute the geodesic distance by doing a breath-first search on the 
Hasse diagram of K(Ti,T2). As we visit each node A in K(Ti,T2), we construct the geodesic 
between T{Xti{A)) and T{A) using the geodesies between T{Xt^{B)) and T{B) for each node B 
covered by A, which we have already visited. As we showed in Section 3, there can be an exponential 
number of elements in the path poset, so this algorithm is exponential in the worst case. However, 
this is a significant improvement over considering each maximal path space. 

Example 5.4. Consider the trees Ti and T2, their incompatibility poset P{Ti,T2), and their path 
poset K{Ti,T2) in Figure |4j Assume the splits in Ti and T2 have lengths, and label each edge in 
the Hasse diagram of K(Ti,T2) with the ratio of the length of the split dropped to the length of 



the split added during the corresponding orthant transition, as shown in Figure 7(a) The sequence 
of ratios along a maximal chain is the ratio sequence passed to PathSpaceGeo to find the path 
space geodesic corresponding to that chain. We first find the geodesic from T({ei,e4}) in O0 to 
T{{fi,f4}) in Oj—j^, by passing the ratio sequences ^,53! and g^f,^ to PathSpaceGeo, 
which returns the distances 1.84 and 1.95, respectively (Figure [7(b) [ ). This and Theorem 5.1 imply 



that if the geodesic from Ti to T2 travels through Oj—^, then it travels through Oj^. We next 
calculate the geodesic between T{{ei, 63, e^}) in O0 and T{{fi, f2, f^}) in Oj—j^. There are three 

maximal chains between and /i, /2. However, we ignore < < /1/4 < /1/2, because we just 
showed that if the geodesic intersects C'j^, then it intersects Oj^. Applying PathSpaceGeo to 
the ratio sequences and gives distances 2.4244 and 2.4243, respectively. 

Thus, by Theorem 5.1 if the geodesic between Ti and T2 intersects Cy^^i then it also intersects 
Oj^, Oj^, and O0. But this is the case, and hence we have found the maximal path space containing 
the geodesic, as shown in Figure [7 (d)[ The length of the geodesic is 2.65. 

We implemented a more memory-efficient version of this algorithm, called GeodeMaps- 
Dynamic. This version uses a depth-first search of the Hasse diagram of K(Ti,T2). For each 
element A in K(Ti,T2), it stores the distance of the shortest path space geodesic found so far 
between T{Xt-^{A)) and T{A). Thus if GeodeMaps-Dynamic revisits an element by following a 
chain in K{Ti,T2) with a longer path space geodesic, it prunes this branch of the search. 

GeodeMaps-Dynamic stores the carrier of the shortest path space geodesic found so far be- 
tween Ti and T2. As a heuristic improvement, at each step in the depth- first search, GeodeMaps- 
Dynamic chooses the node with the lowest transition ratio of the nodes not yet visited. 
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(a) Labelled path (b) First iteration. (c) Second iteration. (d) Third iteration, 
poset K{Ti,T2). 

Figure 7: An example illustrating the dynamic programming approach. 



5.3 GeodeMaps-Divide: a Divide And Conquer Algorithm 

If A is an element in K(Ti,T2), then the trees in Oa share the splits A with T2. This inspires 
the following algorithm, which we call GeodeMaps-Divide. Choose some minimal element of 
P{Ti,T2), and add the splits in this equivalence class to Ti by first dropping the incompatible 
splits. For example, if we choose to add the split set Fi, then we must drop Xtj^(Fi). The trees in 



this orthant Op^ now have splits Fi in common with T2. Apply Theorem 2.1 to divide the problem 
into subproblems along these common splits. For each subproblem, recursively call GeodeMaps- 
Divide. Since some subproblems will be encountered many times, store the geodesies for each 
solved subproblem using a global hash table. 

Example 5.5. Consider the trees Ti and T2 in Figure |8j These trees belong to the family of trees 
given in Figure [5| which have an exponential number of elements in their path posets. Suppose we 
first chose the minimal element /a. We drop 64 from Ti and add /a to get the tree T in Figure 8(c)[ 



T and T2 now share the split fs, so we apply Theorem 2.1 to decompose the problem into two 



subproblems. The incompatibility poset can also be decomposed, as illustrated by Figure 



Each subproblem corresponds to an element in K(Ti,T2), and GeodeMaps-Divide is poly- 
nomial in the number of subproblems solved. Hence an upper bound on the complexity of 
GeodeMaps-Divide is the number of elements in K(Ti,T2). This was shown to be exponen- 
tial in the number of leaves by the family of trees presented in Figure [5} However, for this par- 
ticular family of trees, one can show that GeodeMaps-Divide has a polynomial runtime, while 
GeodeMaps-Dynamic has an exponential runtime. However, there exists a family of trees such 
that GeodeMaps-Dynamic is exponential. See [HI Section 5.2.2] for details. 



5.4 Performance of GeodeMaps-Dynamic and GeodeMaps-Divide 

We now compare the runtime performance of GeodeMaps-Dynamic and GeodeMaps-Divide 
with GeoMeTree [15], the only other geodesic distance algorithm published to our knowledge. For 
n = 10, 15, 20, 25, 30, 35, 40, 45, we generated 200 random rooted trees with n leaves, using a birth- 
death process. Specifically, we ran evolver, part of PAML [30] with the parameters estimated 
for the phylogeny of primates in ^31j, that is 6.7 for the birth rate (A), 2.5 for the death rate {n), 
0.3333 for the sampling rate, and 0.24 for the mutation rate. For each n, we divided the 200 trees 
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Figure 8: An example illustrating the first step of GeodeMaps-Divide. 



into 100 pairs, and computed the geodesic distance between each pair. The average computation 
times are given in Figure [9} Memory was the limiting factor for all three algorithms, and prevented 
us from calculating the missing data points. 

Both GeodeMaps-Dynamic and GeodeMaps-Divide exhibit exponential runtime, but they 
are significantly faster the GeoMeTree. Note that as the trees used were random, they have 
very few common splits. Biologically meaningful trees often have many common splits, resulting 
in much faster runtimes. For example, for a data set of 31 43-leaved trees representing possible 
ancestral histories of bacteria and archaea |16j . we computed the geodesic distance between each 
pair of trees. Using GeodeMaps-Dynamic the average computation time was 0.531 s, while using 
GeodeMaps-Divide the average time was 0.23 s. This contrasts to an average computation time 
of 22 s by GeodeMaps-Dynamic for two random trees with 40 leaves. 

All computations were done on a Dell PowerEdge Quadcore with 4.0 GB memory, and 2.66 
GHz X 4 processing speed. The implementation of these algorithms, GeodeMaps 0.2, is available 
for download from www.cam.cornell.edu/~maowen/geodemaps.html. 

6 Conclusion 

We have used the combinatorics and geometry of the tree space 7^ to develop two algorithms to 
compute the geodesic distance between two trees in this space. In doing so, we have provided a 
linear time algorithm for computing the shortest path in the subspace y(M") of M", which will 
help characterize when the general problem of finding the shortest path through M" with obstacles 
is NP-hard. Furthermore, these algorithms are significantly faster than all known algorithms, and 
the freely available implementation will be of use to any researcher wishing to use the tree space 
framework in their work with phylogenetic trees. For example. Yap and Pachter p2] and Suchard 



23 




Figure 9: Average runtimes of the three geodesic distance algorithms. 



|27j have exphcitly mentioned this as a direction for future work. 



Acknowledgements 

We thank Louis Bihera for numerous helpful discussions and suggestions about this work. We thank 
Karen Vogtmann for sharing her notes and thoughts on the problem, Seth SuUivant for suggestions 
that greatly improved the presentation of this work, and Philippe Lopez for the kind provision of 
the biological data set. Finally, we thank Joe Mitchell for pointing out that finding the geodesic in 
V(R'') is equivalent to solving a touring problem. 



References 

[1] N. Amenta, M. Godwin, N. Postarnakevich, and K. St. John. Approximating geodesic tree 
distance. Inform. Process. Lett., 103:61-65, 2007. 

[2] F. Ardila and C. Klivans. The Bergman complex of a matroid and phylogenetic trees. J. 
Combin. Theory Ser. B, 96:38-49, 2006. 

[3] L. Billera, S. Holmes, and K. Vogtmann. Geometry of the space of phylogenetic trees. Adv. 
in Appl. Math., 27:733-767, 2001. 

[4] G. Birkhoff. Lattice Theory. American Mathematical Society, 1967. 

[5] M.R. Bridson and A. Haefliger. Metric Spaces of Non-positive Curvature. Springer- Verlag, 
1999. 

[6] J. Canny and J. Reif. Lower bounds for shortest path and related problems. In Proceedings 
of the 28th Annual Symposium on Foundations of Computer Science (FOCS), 1987. 

[7] B. DasGupta, X. He, T. Jiang, M. Li, and J. Tromp. On the linear-cost subtree-transfer 
distance between phylogenetic trees. Algorithmica, 25:176-195, 1999. 



24 



[8] B. DasGupta, X. He, T. Jiang, M. Li, J. Tromp, L. Wang, and L. Zhang. Computing distances 
between evolutionary trees. In Handbook of Combinatorial Optimization, pages 35-76. Kluwer 
Academic Publishers, 1998. 

[9] M. Dror, A. Efrat, A. Lubiw, and J. Mitchell. Touring a sequence of polygons. In Proceedings 

of the 35th Annual ACM Symposium on Theory of Computing (STOC), 2003. 

[10] R. Durbin, S. Eddy, A. Krogh, and G. Mitchison. Biological Sequence Analysis: Probabilistic 
Models of Proteins and Nucleic Acids. Cambridge University Press, 1998. 

[11] G.F. Estabrook, F.R. McMorris, and C.A. Meachani. Comparison of undirected phylogenetic 
trees based on subtrees of four evolutionary units. Syst. Zool., 34:193-200, 1985. 

[12] J. Hein. Reconstructing evolution of sequences subject to recombination using parsimony. 

Math. Biosci., 98:185-200, 1990. 

[13] M.D. Hendy and D. Penny. A framework for the quantitative study of evolutionary trees. Syst. 

Zool, 38:297-309, 1989. 

[14] M.K. Kuhner and J. Felsenstein. A simulation comparison of phylogeny algorithms under 
equal and unequal evolutionary rates. Mol. Biol. EvoL, 11:459-468, 1994. 

[15] A. Kupczok, A. von Haeseler, and S. Klaere. An exact algorithm for the geodesic distance 
between phylogenetic trees. J. Comput. Biol, 15:577-591, 2008. 

[16] P. Lopez. Personal communications, 2006. 

[17] J.S.B. Mitchell. Geometric shortest paths and network optimization. In Handbook of Compu- 
tational Geometry, pages 633-701. Elsevier Science, 2000. 

[18] J.S.B. Mitchell and M. Sharir. New results on shortest paths in three dimensions. In 20*'* 
Annual Symposium on Computational Geometry, 2004. 

[19] M. Owen. Distance Computation in the Space of Phylogenetic Trees. PhD thesis, Cornell 
University, 2008. 

[20] V. Polishchuk and J.S.B. Mitchell. Touring convex bodies - a conic programming solution. In 
17th Canadian Conference on Computational Geometry, 2005. 

[21] D.F. Robinson. Comparison of labeled trees with valency three. J. Combinatorial Theory, 
11:105 119, 1971. 

[22] D.F. Robinson and L.R. Foulds. Comparison of phylogenetic trees. Math. Biosci, 53:131-147, 
1981. 

[23] C. Semple and M. Steel. Phylogenetics. Oxford University Press, Oxford, 2003. 

[24] D. Speyer and B. Sturmfels. The tropical Grassmannian. Adv. Geom., 4:389-411, 2004. 

[25] R.P. Stanley. Enumerative Combinatorics, volume 1. Cambridge University Press, 1997. 

[26] A. Staple. Computing distances in tree space. Unpublished research report, Stanford Univer- 
sity, 2004. 

[27] M.A. Suchard. Stochastic models for horizontal gene transfer. Genetics, 170:419-431, 2005. 



25 



[28] K. Vogtmann. Geodesies in the space of trees. Available at 

www.math.cornell.edu/^vogtmann/papers/TreeGeodesicss/index.html, 2007. 

[29] M.S. Waterman and T.F. Smith. On the similarity of dendrograms. J. Theor. Biol., 73:789 - 
800, 1978. 

[30] Z. Yang. PAML 4: a program package for phylogenetic analysis by maximum likelihood. Mol. 
Biol. EvoL, 24:1586-1591, 2007. 

[31] Z. Yang and B. Rannala. Bayesian phylogenetic inference using DNA sequences: A Markov 
Chain Monte Carlo method. Mol. Biol. EvoL, 14:717-724, 1997. 

[32] V.B. Yap and L. Pachter. Identification of evolutionary hotspots in the rodent genomes. 
Genome Research, 14:574-579, 2004. 



26 



